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Abstract. This paper discusses dynamical systems for disk-covering and sphere-packing prob- 
lems. We present facility location functions from geometric optimization and characterize their 
diffcrcntiablc properties. We design and analyze a collection of distributed control laws that are 
related to nonsmooth gradient systems. The resulting dynamical systems promise to be of use in 
coordination problems for networked robots; in this setting the distributed control laws correspond to 
local interactions between the robots. The technical approach relies on concepts from computational 
geometry, nonsmooth analysis, and the dynamical system approach to algorithms. 
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1. Introduction. Consider n sites (pi, . . . ,p„) evolving within a convex polygon 
Q according to one of the following interaction laws: (i) each site moves away from 
the closest other site or polygon boundary, (ii) each site moves toward the furthest 
vertex of its own Voronoi polygon, or (iii) each site moves toward a geometric center 
(circumcenter, incenter, centroid, etc) of its own Voronoi polygon. Recall that the 
Voronoi polygon of the ith site is the closed set of points q g Q closer to pi than to 
any other pj . 

These and related interaction laws give rise to strikingly simple dynamical systems 
whose behavior remains largely unknown. What are the critical points of such dynam- 
ical systems? What is their asymptotic behavior? Are these systems optimizing any 
aggregate function? In what way do these local interactions give rise to distributed 
systems? Does any biological ensemble evolve according to these behaviors and are 
they of any engineering use in coordination problems? These are the questions that 
motivate this paper. 

Coordination in robotics, control, and biology. Coordination problems are 
becoming increasingly important in numerous engineering disciplines. The deploy- 
ment of large groups of autonomous vehicles is rapidly becoming possible because 
of technological advances in computing, networking, and miniaturization of electro- 
mechanical systems. These future multi- vehicle networks will coordinate their actions 
to perform challenging spatially-distributed tasks (e.g., search and recovery opera- 
tions, exploration, surveillance, and environmental monitoring for pollution detection 
and estimation). This future scenario motivates the study of algorithms for autonomy, 
adaptation, and coordination of multi-vehicle networks. It is also important to take 
into careful consideration all constraints on the behavior of the multi- vehicle network. 
Coordination algorithms need to be adaptive and distributed in order for the resulting 
closed-loop network to be scalable, to comply with bandwidth limitations, to tolerate 
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failures, and to adapt to changing environments, topologies and sensing tasks. The in- 
teraction laws introduced above have these properties and, remarkably, they optimize 
network-wide performance measures for meaningful spatially-distributed tasks. 

Coordinated group motions are also a widespread phenomenon in biological sys- 
tems. Some species of fish spend their lives in schools as a defense mechanism against 
predators. Others travel as swarms in order to protect an area that they have claimed 
as their own. Flocks of birds are able to travel in large groups and act as one unit. 
Other animals exhibit remarkable collective behaviors when foraging and selecting 
food. Certain foraging behaviors include individual animals partitioning their envi- 
ronment in nonoverlapping individual zones whereas other species develop overlapping 
team areas. These biological network systems possess extraordinary dynamic capa- 
bilities without apparently following a group leader. Yet, these complex coordinated 
behaviors emerge while each individual has no global knowledge of the network state 
and can only plan its motion according to the observation of its closest neighbors. 

Facility location, nonsmooth stability analysis and cooperative control. 

To analyze the interaction laws introduced above we rely on concepts and methods 
from various disciplines. Facility location problems play a prominent role in the field of 
geometric optimization Facility location pervades a broad spectrum of scientific 

and technological areas, including resource allocation (where to place mailboxes in a 
city or cache servers on the internet), quantization and information theory, mesh 
and grid optimization methods, clustering analysis, data compression, and statistical 
pattern recognition. Smooth multi-center functions for so-called centroidal Voronoi 
configurations and smooth distributed dynamical systems are presented in jX 01 114) . 
Multi-center functions are studied in resource allocation problems ^3 EH] and in 
quantization theory |lf>l 121) . The role of Voronoi tessellations and computational 
geometry in facility location is discussed in [231 HH] ■ 

The notion and computational properties of the generalized gradient are throughly 
studied in nonsmooth analysis [H]- In particular, tools for establishing stability and 
convergence properties of nonsmooth dynamical systems are presented in 1151 H7\ . 
Finally, we refer to I17j for guidelines on how to design dynamical systems for 
optimization purposes, and to 0] for gradient descent flows in distributed computation 
in settings with fixed-communication topologies. 

Recent years have witnessed a large research effort focused on motion planning 
and formation control problems for multi- vehicle systems ^1 1181 1191 1201 1241 1301 ETT] . 
Within the literature on behavior-based robotics, heuristic approaches to the design 
of interaction rules and emerging behaviors have been investigated (see [2] and refer- 
ences therein). Along this specific line of research, no formal results guaranteeing the 
correctness of the proposed algorithms or their optimality with respect to an aggre- 
gate objective are currently available. The aim of this work is to design distributed 
coordination algorithms for dynamic networks as well as to provide formal verifica- 
tions of their asymptotic correctness. A key aspect of our treatment is the inherent 
complexity of studying networks whose communication topology changes along the 
system evolution, as opposed to networks with fixed communication topologies. This 
key aspect is present in the analysis of distributed control laws in O HUE EH and of 
agreement protocols in |24j . 

Statement of contributions. We consider two facility location functions from 
geometric optimization that characterize coverage performance criteria. A collection 
of sites provides optimal service to a domain of interest if (i) it minimizes the largest 
distance from any point in the domain to one of the sites, or (ii) it maximizes the 



Coordination and geometric optimization via distributed dynamical systems 



3 



minimum distance between any two sites. In other words, if P = (pi, . . . ,p n ) are n 
sites evolving within a convex polygon Q, we extremize the multi-center functions 

max min d(q,Pi) > , min {h d {Pi,Pj), d(p i: dQ)\ , 

q€Q |_ie{l,...,n} J 

where d(p, q) and d(p, <9Q) are the distances between p and g, and between p and 
the boundary of Q, respectively. (The role of the \ factor will become clear later.) 
We study the differentiable properties of these functions via nonsmooth analysis. We 
show the functions are globally Lipschitz and regular, we compute their generalized 
gradients, and we characterize their critical points. Under certain technical condi- 
tions, we show that the local minima of the first multi-center function are so-called 
circumccnter Voronoi configurations, and that these critical points correspond to the 
solutions of disk-covering problems. Similarly, under analogous technical conditions, 
we show that the local maxima of the second multi-center function are so-called incen- 
ter Voronoi configurations, and that these critical points correspond to the solutions 
of sphere-packing problems. 

Next, we aim to design distributed algorithms that extremize the multi-center 
functions. Roughly speaking, by distributed we mean that the evolution of each site 
depends at most on the location of its own Voronoi neighbors. We study the gener- 
alized gradient flows induced by the multi-center functions using nonsmooth stability 
analysis. Although these dynamical systems possess some convergence properties, 
they are not amenable to distributed implementations. Next, drawing connections 
with quantization theory, we consider two dynamical systems associated to each multi- 
center function. First, we consider a novel strategy based on the generalized gradient 
of the 1-ccnter functions of each site, and, second, we consider a geometric centering 
strategy similar to the well-known Lloyd algorithm I^T] . 

Remarkably, these strategies arising from the nonsmooth gradient information 
have natural geometric interpretations and are indeed the local interaction rule de- 
scribed earlier. For the first (respectively second) multi-center function, the first 
strategy corresponds to the interaction law "move toward the furthest vertex of own 
Voronoi polygon" (respectively, "move away from the closest other site or polygon 
boundary", and the second strategy corresponds to the interaction law "move to- 
ward circumcenter of own Voronoi polygon" (respectively "move toward incenter of 
own Voronoi polygon"). We prove the uniqueness of the solutions of the resulting 
distributed dynamical systems and we analyze their asymptotic behavior using nons- 
mooth stability analysis, showing that the active sites will approach the corresponding 
centers of their own Voronoi cells. 

Two of our results are related to well-known conjectures in the locational opti- 
mization literature |13l I29| : (i) that the first multi-center problem is equivalent to 
a disk covering problem (how to cover a region with possibly overlapping disks of 
equal minimum radius), and (ii) that the generalized Lloyd strategy "move toward 
circumcenter of own Voronoi polygon" converges to the set of circumcenter Voronoi 
configurations. 

Organization. The paper is organized as follows. Section [5] provides the pre- 
liminary concepts on Voronoi partitions, nonsmooth analysis, stability analysis, and 
gradient flows, and introduces the multi-center problems. Section presents a com- 
plete treatment on the functions analysis and algorithm design for the 1-center prob- 
lems. Section 0] discusses the differentiable properties and the critical points of the 
multi-center functions. Section [S] introduces a number of dynamical systems (smooth 
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and nonsmooth, distributed and non-distributed) and analyzes their asymptotic cor- 
rectness. 

2. Preliminaries and problem setup. Let || • j| denote the Euclidean distance 
function on and let v • w denote the scalar product of the vectors v, w 6 R w . Let 
vrs(ii) denote the unit vector in the direction of ^ « £ R , i.e., vrs(w) = w/IMI- 
Given a set S in l", we denote its convex hull by co^) and its interior set by 
int(,S). If S is a convex set in 1^, let proj s : R N — > S denote the orthogonal 
projection onto S and let D5 : R N — > R denote the distance function to S. For 
R > 0, let B N (p,R) = {qeR N \ \\p-q\\ < R}, and B N (p,R) = mt(S N (p, R)). A 
set {vi, . . . , vm} of vectors in R N positively spans R N if any w G M. N can be written 
as w = 'YldLi a l v h with ai > 0, I £ {1, . . . , M}. The following simple lemma, e.g., 
see [S], characterizes this situation. 

Lemma 2.1. Given a set {vi, . . . ,vm} of AI arbitrary vectors in R N , then the 
following statements are equivalent 

(i) ...,%} positively spans R N ; 

(ii) G int(co{wi, . . .,%}); 

(Hi) for each w G R N , there exists Vi such that w ■ Vi > 0. 

Let Q be a convex polygon in R 2 . We denote by Ed(Q) = {ei, . . . , ejvf} and 
Ve(Q) = {vi, . . . ,vl} the set of edges and vertexes of Q, respectively. Let P = 
(pi, . . . ,p n ) G Q n C (R 2 ) n denote the location of n generators in the space Q. Let 
TTi : Q" —> Q be the canonical projection onto the ith factor, 7Tj(pi, . . . ,p n ) = pi. Note 
that this mapping is surjective, continuous and open (the latter meaning that open 
sets of Q n arc mapped onto open sets of Q). 

2.1. Voronoi partitions. We present here some relevant concepts on Voronoi 
diagrams and refer the reader to |11II23| for comprehensive treatments. A partition of 
Q is a collection of n polygons W = {Wx, . . . , W n } with disjoint interiors whose union 
is Q. Of course, more general types of partitions could be considered (as, for instance, 
continuous deformations of the previous ones) , but these ones will be sufficient for our 
purposes. The Voronoi partition V(P) = (Vi(P), . . . ,V n (P)) of Q generated by the 
points (pi, . . . ,p n ) is defined by: 

Vi(P) = {qeQ\ \\q-Pi\\ < ||? -ft || , Vj ? 1} . 

For simplicity, we shall refer to Vi(P) as Vi. Since Q is a convex polygon, the boundary 
of each Vi is the union of a finite number of segments. If Vi and Vj share an edge, 
i.e., V, H V} is neither empty nor a singleton, then pi is called a (Voronoi) neighbor 
of pj (and vice- versa) . All Voronoi neighboring relations are encoded in the mapping 
Af : Q n X {1, . . . ,n} — > 2* 1 ' ' ,n ^ where M(P,i) is the set of indexes of the Voronoi 
neighbors of p;. Of course, j G J\f(P,i) if and only if i £ J\f(P,j). We will often omit 
P and instead write M(i). 

For P G Q n , the vertexes of the Voronoi partition V(P) are classified as follows: 
the vertex v is of type (a) if it is the center of the circle passing through three gener- 
ators (say, pi, pj, and pk), the vertex v is of type (b) if it is the intersection between 
an edge of Q and the bisector determined by two generators (say, e, pi, and pj), and 
the vertex v is of type (c) if it is a vertex of Q, i.e., it is determined by two edges of 
Q and by the generator of a cell containing it (say, e, /, and pi). Correspondingly, we 
shall write v(i,j, k), v(e, and v(e, f, i) respectively, whenever we are interested in 
making explicit the elements defining the vertex v. The vertex v G Ve(V,(P)) is said 
to be nondegenerate if it is determined by exactly three elements (e.g., as described 
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above, either three generators, or an edge and two generators, or two edges and one 
generator), otherwise it is said to be degenerate. Further, the configuration P is said 
to be nondegenerate at the ith generator if all vertexes v €E Ve(Vi(P)) are nondcgcn- 
erate, otherwise P is degenerate at the ith generator. Finally, a configuration P is 
said to be nondegenerate if all its vertexes are nondegenerate, otherwise it is said to 
be degenerate. These concepts are illustrated in Fig. 12.11 
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Fig. 2.1. A Voronoi partition with degenerate and nondegenerate vertexes. Vertexes v a , t>6, and 
v c are nondegenerate vertexes of type (a), (b), (c), respectively. Vertexes Vd and v e are degenerate. 

For P € Q n , the edges of the Voronoi partition V(P) are classified as follows: the 
edge e is of type (a) if it is a segment of the bisector determined by two generators 
(say, Pi, pj), the edge e is of type (b) if it is contained in the boundary of Q, i.e., 
it is a subset of an edge of Q and it belongs to a single cell (say, the cell of the 
generator j»i). Correspondingly, we shall write e(i,j) and e(i) respectively, whenever 
we arc interested in making explicit the elements defining the edge e. Further, when 
considering an edge of type (a), we let n e uj) denote the unit normal to e(i,j) pointing 
toward int(Vi(P)). When considering an edge of type (b), we let n e /j) denote the unit 
normal to e(i) pointing toward int(Q). 

2.2. The disk-covering and the sphere-packing problems. We are inter- 
ested in the following locational optimization problems 



min max min \\q - Pi\\ } } , (2.1) 

Pi,—,Pn I gGQ I ie{l,...,n} 



max < min {h\\Pi ~ Pj\\,^>e(pi)} 

Pl,—,Pn I ij"G{l,...,n} 

ijtj, eGEd(Q) 



(2.2) 



The optimization problem (|2.1|) is referred to as the p-center problem in |131 129j . 
Throughout the paper, we will refer to it as the multi-circumcenter problem. In the 
context of coverage control of mobile sensor networks |1(J| . the multi-circumcenter 
problem corresponds to considering the worst case scenario, in which no information 
is available on the distribution of the events taking place in the environment Q. The 
network therefore tries to minimize the largest possible distance of any point in Q to 
one of the generators' locations given by p±, . . . ,p n , i.e. to minimize the function, 



TLdc (P) — max s m hi \\l — Pi\\f= max < max 1 1 o — ». 
qeQ {ie{i,...,n} J ie{l,...,n} {qeVi 

It is conjectured in |29j that this problem can be restated as a disk-covering problem: 
how to cover a region with (possibly overlapping) disks of minimum radius. The 
disk-covering problem then reads: 

min{i? | U l6 {i i ... i „} B 2 {pi,R) 2 Q} ■ 
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Wc shall present a proof of this statement in Theorem 14.71 below. Given a polytopc 
W in M. N , its circumcenter, denoted by CC(W), is the center of the minimum-radius 
sphere that contains W. The circumradius of W, denoted by CR(W), is the radius 
of this sphere. Wc will say that P is a circumcenter Voronoi configuration if pi = 
CC(Vi(P)), for all i G {1, . . ., n}. We denote by Vc DC (V(P)) the set of vertexes of the 
Voronoi partition where the value TCbc(P) is attained, i.e. v G VeDc(V(P)) if there 
exists i such that v G Vi{P) and \\v — Pi\\ = Hdc(P)- 

We will refer to the optimization problem (|2.2|) as the multi-incenter problem. In 
the context of applications, this problem corresponds to the situation where we are 
interested in maximizing the coverage of the area Q in such a way that the sensing 
radius of the generators do not overlap (in order not to interfere with each other) or 
leave the environment. We therefore consider the maximization of the function 

T-isp(P) = . min -pjUDeipi)} = min <^ min \\q-pi 

iS{l,...,n} ^qgint(Vi) 

i#j,e6Ed(Q) 

A similar conjecture to the one presented above is that the multi-incenter problem 
can be restated as a sphere-packing problem: how to maximize the coverage of a 
region with non-overlapping disks (contained in the region) of minimum radius. The 
problem reads: 

max{P| U, 6(1 „)B 2 (p„JJ)C(3, B 2 ( Pi , R) n B 2 (p j , R) = 0} . 

In Theorem 14.81 we provide a positive answer to this question. Given a polytope W 
in R N , its incenter set (or Chebyshev center set, see |BJ), denoted by IC(W), is the 
set of the centers of maximum-radius spheres contained in W. The inradius of W, 
denoted by IR(VT'), is the common radius of these spheres. We will say that P G Q n 
is an incenter Voronoi configuration if pi G lC(Vi(P)), for all i G {1, . . . , n}. If P is an 
incenter Voronoi configuration, and each Voronoi region Vi(P) has a unique incenter, 
lC(Vi(P)) = {pi}-, then wc will say that P is a generic incenter Voronoi configuration. 
We denote by Edsp(V(P)) the set of edges of the Voronoi partition where the value 
Wsp(-P) is attained, i.e. e G Edsp(V(P)) if there exists i such that e G Ed(Vi(P)) and 

D e fe) =«SP(P). 

2.3. Nonsmooth analysis. The following facts on nonsmooth analysis |Hj will 
be most helpful in analyzing the properties of the locational optimization functions 
for the disk-covering and the sphere-packing problems, as well as the convergence of 
the distributed algorithms wc will propose to extremize them. 

Definition 2.2. A function f : R N — > R is said to be locally Lipschitz near 
x G M. N if there exist positive constants L x and e such that \f(y) — f(y')\ < Lx\\y — y'\\ 
for all y,y' G B N (x,e). 

Note that continuously differentiable functions at x are locally Lipschitz near x. 
The usual right directional derivative of / at x in the direction of v G R N is defined 

MS 

v r f(x + tv)-f(x) 

f (x,V) = hm , 

t^o+ t 

when this limits exists. On the other hand, the generalized directional derivative of / 
at x in the direction of v G R N is defined as 

j (x; v) = hm sup . 

y— >x t 
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This notion of directional derivative has the advantage of always being well-defined. 

Definition 2.3. A function f : R N — > R is said to be regular at x £ M. N if for 
all v £ M. N , f'{x; v) exists and f°(x; v) = f'(x; v). 

Again, a continuously diffcrcntiablc function at x is regular at x. Also, a locally 
Lipschitz function at x which is convex (or concave) is also regular (cf. Proposi- 
tion 2.3.6 in EH). 

From Radcmachcr's Theorem we know that locally Lipschitz functions are 
continuously differentiable almost everywhere (in the sense of Lebesgue measure) . If 
O f denotes the set of points in R w at which / fails to be differentiable, and S denotes 
any other set of measure zero, the generalized gradient of / is defined by 



Note that this definition coincides with df(x) if / is continuously differentiable at 
x. The generalized gradient and the generalized directional derivative (cf. Proposi- 
tion 2.1.2 in [§]) are related by f°(x;v) = max{C • v \ ( £ df(x)}, for each v £ R N . 
A point x £ R w which verifies that £ df(x) is called a critical point of f. 
The following result corresponds to Proposition 2.3.12 in [§]. 

PROPOSITION 2.4. Let {f k : R N — > R | k £ {1, . . . , m}} be a finite collection of lo- 
cally Lipschitz functions nearx £ H N . Consider f(x') — min {fk(x') \ k E {1, . . . , to}}. 
Then, 

(i) f is locally Lipschitz nearx, 

(ii) if I(x') denotes the set of indexes k for which fk(x') = f(x'), we have, 



and if each fi is regular at x for i £ I(x), then equality holds and f is regular 
at x. 

The extrema of Lipschitz functions are characterized by the following result. 

Proposition 2.5. Let f be a locally Lipschitz function at x £ R N . If f attains 
a local minimum or maximum at x, then £ df(x), i.e., x is a critical point. 

Let Ln : 2 R ™ — > 2 R ™ be the set- valued mapping that associates to each subset S 
of H N the set of its least- norm elements Ln(S). If the set S is convex, then the set 
Ln(5) reduces to a singleton and we note the equivalence Ln(5) = proj s (0). Along 
the paper, we shall only apply this function to convex sets. For a locally Lipschitz 
function /, we consider the generalized gradient vector field Ln(0/) : M. N — > M. N given 



Theorem 2.6. Let f be a locally Lipschitz function at x. Assume df(x). 
Then, there exists T > such that 



The vector — Ln(df)(x) is called a direction of descent. 

2.4. Stability analysis via nonsmooth Lyapunov functions. For differen- 
tial equations with discontinuous right-hand sides, solutions are defined in terms of 
differential inclusions |15| . 

Let F : R N -> 2 R be a set- valued map. Consider the differential inclusion 




df(x) C co {dfi(x) \i£l(x)} , 



(2.3) 



by x i— > 



Ln(0/)(a:) = Ln(0/(a:)). 



/(x-iLn(9/)(a;))</(x)-^|Ln(9/)(a:)|| 2 , 0<<<T. 



x £ F(x) 



(2.4) 



8 



Jorge Cortes and Francesco Bullo 



A solution to this equation on an interval [to,ii] C R is defined as an absolutely 
continuous function x : [io,ti] — * R N such that x(t) G F(x(t)) for almost all t € 
[io,ti]. Given xq G R^, the existence of at least a solution with initial condition xn 
is guaranteed by the following lemma. 

Lemma 2.7. Let the mapping F be upper semicontinuous with nonempty, compact 
and convex values. Then, given xq G M. N , there exists at least a solution of 1)2.4(1 with 
initial condition Xq- 

Now, consider the differential equation 

x(t)=X(x(t)), (2.5) 

where X : M. N — > M. N is measurable and essentially locally bounded. The solution of 
this equation has to be understood in the Filippov sense. For each x G M. N , consider 
the set 

K[X](x)=f] p| co{X(B N (x,S)\S)}, 

5>0^(S)=0 

where \x denotes the usual Lebesgue measure in M. N . Alternatively, one can show [2U 
that there exists a set Sx of measure zero such that 

K[X](x) = co\ lim X(xi) \ £ S U S x > , 

>+oo J 

where S is any set of measure zero. A Filippov solution of l|2.5|l on an interval 
[to, ti] C K is defined as a solution of the differential inclusion 

xGK[X](x). (2.6) 

Since the multivalued mapping : M N —> 2 R is upper semicontinuous with 

nonempty, compact, convex values and locally bounded (cf. |15|). the existence of 
Filippov solutions of (|2.5fl is guaranteed by Lemma \l. 71 

A set M is weakly invariant (respectively strongly invariant) for l|2.5|l if for each 
xq G M, M contains a maximal solution (respectively all maximal solutions) of l|2.5|) . 
Given a locally Lipschitz function / : M. N — > R, the set-valued Lie derivative of f with 
respect to X sX x is defined as 

Cxf{x) = {a G R | 3v G K[X]{x) such that Q ■ v = a , V£ G df{x)} . 

For each x G M. N , Cxf(x) is a closed and bounded interval in R, possibly empty. 
If / is continuously differentiable at x, then Cxf(x) = {df ■ v | v G A'[A](x)}. If, in 
addition, X is continuous at x, then Cxfix) corresponds to the singleton {Cx f{x)}, 
the usual Lie derivative of / in the direction of X at x. The importance of the 
set- valued Lie derivative stems from the next result 

Theorem 2.8. Let x : [io,<i] — > R w be a Filippov solution of l|2.5|l . Le£ f be a 
locally Lipschitz and regular function. Then 4r (f(x(t))) exists a.e. and -i (f(x(t)j) G 

C x f(x(t)) a.e. 

The following result is a generalization of LaSalle principle for differential equa- 
tions of the form 1(2.5(1 with nonsmooth Lyapunov functions. The formulation is taken 
from [3], and slightly generalizes the one presented in |27( . 

Theorem 2.9 (LaSalle principle). Let f : R N — > R be a locally Lipschitz and 
regular function. Let xq G M. n and let / _1 (< f(xo),Xo) be the connected component of 
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{x G M. N | f(x) < /(ceo)} containing xq. Assume the set f 1 (< f(xo),Xo) is bounded 

and assume either maxCx f(x) < or Cxfix) = for all x G / _1 (< f(xo),Xo)- 
Then / (< f(xo),Xo) is strongly invariant for <|2.5[) . Let 

Z x ,f = G R N | G Cxf(x)} . 

Then, any solution x : [to,+oo] — > M. N of <|2.5[) starting from xo converges to the 
largest weakly invariant set M contained in Zxj H / _1 (< f(xo),Xo). Furthermore, 
if the set M is a finite collection of points, then the limit of all solutions starting at 
xq exists and equals one of them. 

The proof of the last fact in the theorem statement is the same as in the smooth 
case, since it only relies on the continuity of the trajectory. The next statement is 
based on Theorem 2 of |25| . 

Proposition 2.10. Under the same assumptions of Theorem \2.iA !/max£x/(i) < 
— e < a.e. on R \ Zxj, then Zxj is attained in finite time. 

Proof. Let x : [to, +oo) — > R N be a Filippov solution starting from xq- We argue 
that there must exist T such that x(T) G Zx.f- Otherwise, we have 

f(x(t)) = f(x(t )) + f ±f(x(s))ds < f(x(t )) - e{t - t ) -co , 

Jto ds 

which contradicts the fact that /~ 1 (< f(xo),Xo) is strongly invariant and bounded. 
□ 

2.5. Nonsmooth gradient flows. Finally, we are in a position to present the 
nonsmooth analogue of well-known results on gradient flows. Given a locally Lipschitz 
and regular function /, consider the following generalized gradient flow 

x(t) = -Ln(df)(x(t)). (2.7) 

Theorem 12.61 guarantees that, unless the flow is at a critical point, —Ln(df)(x) is 
always a direction of descent at x. In general, the vector field Ln(df) in l|2.7|l is 
discontinuous, and therefore its solution must be understood in the Filippov sense. 
Note that, since / is locally Lipschitz, Ln(df) = df almost everywhere. An important 
observation in this setting is that K[df](x) = df(x) (cf. [21] )■ The following result, 
which is a generalization of the discussion in jS] , guarantees the convergence of this 
flow to the set of critical points of /. 

PROPOSITION 2.11. Let x G R n and assume / _1 (< f(x ),x ) is bounded. 
Then, any solution x : [to, +oo) — ► R N of eq. I|2.7|) starting from xq converges asymp- 
totically to the set of critical points of f contained in /~ 1 (< f(xo),Xo). 

Proof. Let a G £_L n (a/)/(a;)- By definition, there exists w G K[— Ln(df)](x) = 
—df(x) such that a = w ■ ( for all £ G df(x). In particular, for ( = —w G df(x), 
we have a = — \\w\\ 2 < 0. Therefore, max£_ Ln (a/)/(i) < or £_ hn(df)f( x ) = 
Now, resorting to the LaSalle principle fTheorem l2.9|l . we deduce that any solution 
x : [to,+oo) — * R N starting from xo converges to the largest weakly invariant set 
contained in •^_Ln(o/),/ H / _1 (< f(xo),xo). Let us see that -^_Ln(a/),/ is equal 
to L = {x G Q n G df(x)}. Obviously, L C ^-Ln(9/),/- Conversely, assume 

x G ^-Ln(9/),/- Then, G Ln(df)f( x ), i- e -i there exists v G —df(x) such that 
£ • v = for all C G df(x). In particular, for ( = —v, we get ||w|| 2 = 0, that 
is, v = G df(x), as desired. Note that ^-Ln(a/),/ = Lq is the equilibrium set 
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of i|2.7[) and therefore is weakly invariant. Finally, we prove that it is also closed. 
Let x £ ^-Ln(a/),/ an d consider a sequence {x^ £ R N | k £ N} C ^-Ln(a/)./ such. 
that Xk — > x. Then, using the fact that the multivalued mapping K[— v] is upper 
scmicontinuous, for any e > 0, there exists fco such that for k > ko, df(xk) C df(x) + 
Sjv(0, e). Since Xk £ 2_ Ln (a/)j, then £ df(x) + Bn(0,e) for all e > 0, and this 
implies that E df(x), i.e., x E ^-Ln(a/)./- Hence the largest weakly invariant 
set contained in Z_L n (e/),/ n / _1 (< /(zo),£o) is ^-Ln(a/),/ H / _1 (< /(^o),^o) = 
{^e/- 1 ^/^),^) |0ES/(x)}. □ 

3. The 1-center problems. In this section we consider the disk-covering and 
the sphere-packing problems with a single generator, i.e., n= 1. This treatment will 
give us the necessary insight to tackle later the more involved multi-center version of 
both problems. When n = 1, the minimization of Tine simply consists of finding the 
center of the minimum-radius sphere enclosing the polygon Q. On the other hand, the 
maximization of Tisp corresponds to determining the center of the maximum-radius 
sphere contained in Q. Let us therefore define the functions 

\Eq(p) = max{||?-p|| \ q£Q} = max{||i> -p|| | v £ Ve(Q)} , 
sm Q (p) =min{||?-p|| | q & int(Q)} = min{D e (p) | e E Ed(Q)} . (3.1) 

When n = 1, we then have that Hdc = IgQ : Q - * R an d 7"fep = sm Q : Q ^ M. 

3.1. Smoothness and critical points. We here discuss the smoothness prop- 
erties and the critical points of the 1-center functions. Since the function lgg is the 
maximum of a (finite) set of convex functions in p, it is also a convex function 0. 
Therefore, any local minimum of lgg is also global. 

Lemma 3.1. The function lgg has a unique global minimum, which is the cir- 
cumcenter of the polygon Q. 

Proof. Let F : R — > R be any continuous non-decreasing function. Then, 

F(lg Q (p))=max{F(\\v-p\\)\v£Ve(Q)} . 

If we take F(x) = x 2 , each function \\v — p\\ 2 is strictly convex, and hence F(lgg(p)) 
is also strictly convex. Therefore, this latter function has a single minimum on Q. 
Since any global minimum of lgg is also a global minimum of F(lgQ(p)), we conclude 
the result. □ 

The function sniQ is the minimum of a (finite) set of affine (hence, concave) 
functions defined on the half-planes determined by the edges of Q, and hence it is 
also a concave function [S] on the intersection of their domains, which is precisely Q. 
Therefore, any local maximum of smQ is also global. However, this maximum is not 
unique in general. 

Lemma 3.2. The incenter set of the polygon Q is the set of maxima of the 
function sitiq and it is a segment. 

Proof. It is clear that the set of maxima of smQ is IC(Q). As a consequence of 
the concavity of smQ over the convex domain Q, one deduces that IC(Q) is a convex 
set. Now, assume there are three points Pi,P2,P3 in IC(Q) which are not aligned. 
Since B 2 (q,lR(Q)) C Q for all q £ co(pi,p 2 ,P3) C IC(Q), and co(pi,p 2 ,P3) has non- 
empty interior, there exist qo £ Q and r > IR(<5) such that B2{qo, r) C Q, which is a 
contradiction. □ 

Note that the circumcenter of a polygon can be computed via the finite-step 
algorithm described in |28j . The incenter set of a polygon can be computed via the 
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following linear program: maximize the inradius subject to the constraints that the 
distance between the incenter and each of the polygon edges must be greater than 
or equal to the inradius. In what follows, let us examine dynamical systems that 
compute these geometric centers. 

Proposition 3.3. The functions Igg(p), smg(p) are locally Lipschitz and regu- 
lar, and their generalized gradients are given by 

d\g Q {p) = co{vrs{p-v)\veVe{Q), \g Q {p) = \\p - v\\} , (3.2) 

dsm Q (p) = co{n e \ e £ Ed(Q) , sm Q (p) = D e (p)} • (3.3) 

Moreover, 

e dlg Q ( P ) ^ p = CG(Q) , Oedsm Q (p)^pelC(Q), (3.4) 

and, i/0 S int(<9smg(p)), then IC(Q) = {p}. 

Proof. Given the expressions in l|3.1|l and Proposition 12.41 we deduce that lgg 
and smg are locally Lipschitz and regular, and that their generalized gradients are 
respectively given by (|3.2|l and l|3.3|l . Concerning (|3.4|) . the implications from right to 
left in H3.4fl readily follow from Proposition ^. 51 As for the other ones, note that it is 
sufficient to prove that p is a local minimum, respectively that p is a local maximum. 
We prove the result for the function lgg. The proof for smg is analogous. Assume that 

G d\gq(p). Then, there exist vertexes , . . . , Vi K of Q with lgg(p) = \\vi t — p\\, 

1 G {1,...,K} such that = J2ie{i....,K} ^ VTS (P ~ w »i)> where J2ie{i,...,K} ^ = L 
A; > 0, I G {1, . . . , K}. Let U be a neighborhood of p and take q £ U . One can 
show that there must exist I* such that (p — v^, ) ■ (q — p) > 0, since otherwise 
= • (q — p) = (X^efi k} ^ l VTS (p ~ v ii)) ' (q ~ p) < 0, which is a contradiction. 
Then, 

Ik - f = h -P\\ 2 + \\P- V it . f - 2(q - p) ■ (v h , -p)> \\p - v H , f . 

Therefore, lgg(g) > ||p — v^, \\ = lgg(p), which shows that p is a local minimum. 
Finally, if G int(9 snig(p)), then one can see that p is a strict local maximum. 
Furthermore, there cannot be any other local (hence global) maximum of smg, as we 
now show: assume p G IC(<5). By hypothesis, the sphere B 2 (p, smg(p)) centered at p 
of radius smg(p) is contained in Q. Consider the vector p — p. By Lemma 12. II there 
exists e G Ed(Q) with T> e (p) = smg(p) such that (p—p)-n e > 0. Therefore, there are 
points of Bi(j>, smg(p)) which necessarily belong to the half-plane defined by e where 
Q is not contained, which is a contradiction. □ 

3.2. Convergence properties for nonsmooth gradient flows. Here we study 
the generalized gradient flows arising from the two 1-center functions. An immediate 
consequence of Propositions 12. 1 fl and 13.31 is the following result. 

Corollary 3.4. The gradient flows of the functions lgg and smg 

x(t)=-Ln(dlg Q )(x(t)), (3.5) 
x(t) = Ln(d smg) (#(*)) , (3.6) 

converge asymptotically to the circumcenter CC(Q) and the incenter set IC(Q), re- 
spectively. 

The following two propositions discuss the convergence properties of the gradient 
descents. 

Proposition 3.5. IfO e int(<91gg(CC(Q))), then the flow (££3 reaches CC(Q) 
in finite time. 
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Proof. Let us prove that there exists e > such that max£_ Ln [ lgQ ] lgg < — e a.e. 
on Q \ {CC(Q)}. Take p ^ CC(Q). We know that each element a G £-Ln[i gQ ] ^Eq(p) 
can be expressed as a = — ||u>|| 2 , with —we d\gq{p). Therefore, we have 

maxL Ln[lgQ] lg Q (p) = -|| Ln[lg Q ](p)|| 2 . 

If there is a single vertex of Q involved in d\gq{p), then moving along the direction 
— Ln[lgQ](p) obviously decreases the distance to that vertex while maintaining con- 
stant the norm of the least-norm element, which is 1. If there are two or more vertexes 
involved, from the expression for the generalized gradient at p (cf. eq. I|3.2|l '). it is 
clear that one can express it as 

{xeR N \gi(x)<0,..., 9s (x)<0} , 

for some linear functions g r . Note that the points x € <91gg(p) such that g r {x) = 
for some r correspond to a set of the form co{vrs(p — w r .i), vrs(p — iv.2)}, for certain 
vertexes v r ,i, tV,2 of Q. Now, the computation of the least-norm element in d\gg(p) 
can be formulated as the convex problem, 

minimize ||a;|| 2 

subject to gi(x) < 0, . . . , g s {x) < . 

Let x* = Ln[lgg](p). Let R denote the set of indexes r for which g r (x*) = 0. Then x* 
is a regular point [22], meaning to say that dg r (x*), r S R are linearly independent 
vectors. This is because the cardinality of R is at most 2 (since the intersection 
of two lines already determines a point), and the gradients of the functions g r are 
independent when considered pairwise. We apply then the Kuhn- Tucker first-order 
necessary conditions for optimality [22] to conclude that there must exist r* £ R such 
that g r *(x*) = 0. It is easy to sec that r* must be unique, since otherwise x* does 
not have minimum norm. Therefore, we have that Ln[lgg](p) is determined as the 
least-norm element in co{vrs(p — iv*.i),vrs(p — W r *,2)}- As a consequence, moving 
along the direction — Ln[lgg](p) decreases the distance to the vertexes tv*,i, v r * t 2, 
and hence the norm of the least-norm element decreases. If, along the flow l|3.5f) . a 
new vertex of Q enters in the computation of <91gg(p(i)), then there can be a jump in 
the norm of Ln[lgQ](p(£)), which by definition will always be decreasing. Finally, note 
that if iv*, ij IV*, 2 are active at the circumcenter, then they cannot be opposite with 
respect to CC(Q) precisely because of the assumption that lies in int (<9 lgg(CC(Q))). 
Therefore, we conclude 

|| Ln[lg Q ](p)|| > e = min {1, {|| Ln(co{vrs(CC(Q) - v), vrs(CC(Q) - w)})\\ \ 

v,we 7(CC(Q)), CC(Q) -v? -(CC(Q) - w)}} > , Vp ^ CC(Q) . 

Resorting now to Proposition ^. 101 we deduce that the circumcenter CC(Q) is attained 
in finite time. □ 

Remark 3.6. Note that if G dlg Q (CC(Q))\mt(<91g Q (CC(Q))), then generically 
convergence is achieved over an infinite time horizon. 

Proposition 3.7. The flow l|3.6|l reaches the set IC(Q) in finite time. 

Proof. Let p ^ IC((J). We know min£ Ln [ smQ ] smg(p) = || Ln[sniQ] (p) || 2 . More- 
over, for all p g" IC(<5), we have 

||Ln[sm Q ](p)|| > e = min{l,{||Ln(co{n e ,n / })|| | e, / £ Ed(Q),n e ^ -n f }}> 0. 
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Resorting to Proposition ^. 1UI wc deduce the desired result. □ 

Fig. 13.11 shows an example of the implementation of the gradient descent (I3.5|l 
and (|3.6() . Note that if the circumcenter CC(Q) (respectively the incenter set IC(Q)) is 
first computed offline, then the strategy of directly going toward it would converge in a 
less "erratic" way. Note also that the move-toward-the-center strategy is exponentially 
fast. 




Fig. 3.1. Illustration of the gradient descent of lgg and smQ. The points where the curve 
t \— * p(t) fails to be differentiable correspond to points where there is a new vertex v of Q such 
that \\p(t) — v\\ = lgQ(p(t)) (respectively a new edge e of Q such that D e (p(t)) = smQ(p(t))). The 
circumcenter and the incenter are attained infinite time according to Provositions VJ .E\ and X$7\ 

Finally, we conclude this section with four results useful for later developments. 
Lemma 3.8. Let q £ Q, let v(q) be one of the vertexes of Q which is furthest 
away from q, and let e(q) be one of the edges of Q which is nearest to q. Then, 
(i) Ln[lgg](g) ■ {q — v(q)) > 0, and the inequality is strict if q ^ CC(Q), 
(ii) (q - CC(Q)) • (q - v(q)) > \\q - CC(Q)|| 2 /2, 

(Hi) Ln[smg](q) • n e > 0, and the inequality is strict if q $ IC(<5), and 

(iv) (x — q)-n e > IR(Q) — D e (g) > for any x E IC(Q), and the second inequality 

is strict if q £" IC(Q). 
Proof. Let q be a point in Q. If q = CC(Q), claims (i) and (ii) arc obviously 
satisfied since Ln[lgg](g) = 0. Assume then that q ^ CC(Q). Let us prove first 
(i). By Proposition 13.31 ^ <91gg(q), and hence Ln[lgg](g) ^ 0. Let us prove 
Ln[lgg](q) • (q — v(q)) > reasoning by contradiction. If Ln[lgg](g) • (q — v(q)) < 0, 
then d/dt(\\q-tLn[lg Q ](q) - v\\) t=Q = vrs(q - v) ■ (- Ln[lg Q ](g)) > 0, which im- 
plies that \\q — tLn[lgg](g) — v\\ > \\q — v\\ = lgg(<;) for t > small enough. On 
the other hand, invoking Theorem 12.61 we have that lgg(g) — t\\ Ln[lgg](g)|| 2 /2 > 
1§q(9 — t Ln[lgg](g)) > \\q — tLn[lgg](q) — u||. Gathering both facts, we conclude 
—t\\ Ln[lgg](q)|| /2 > 0, which is a contradiction. 

Let us now prove (ii). Since q ^ CC(<3), we have \\q — v(q)\\ > CR(Q). Consider 
then a ball B2{v(q), \\q — w(<?)||) centered at the vertex v(q), with radius \\q — v(q)\\. 
By definition of the circumcenter, CC(Q) must lie in the interior of B2(v(q), \\q — 
v(q)\\). Consequently, || CC(Q) - v(q)\\ < \\q-v(q)\\. Then, from || CC(Q) - v(q)\\ 2 = 
|| CC(Q) -q\\ 2 + \\q - v(q)\\ 2 - 2(q - CC(Q)) • (q - v(q)), we deduce 

2(q - CC(Q)) • (q - v(q)) - \\ CC(Q) - q\\ 2 = \\q - v(q)\\ 2 - \\ CC(Q) - v(q)\\ 2 > , 



which implies the desired result. 
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Let us now prove (iii). If g G IC(Q), the claim is obviously satisfied since 
Ln[smg](g) = 0. Assume then that q IC(Q). By Proposition ^. 31 £" dsuiQ(q), and 
hence Ln[sm.Q](g) ^ 0. Let us prove Ln[smg](q) ■ n e > reasoning by contradiction. 
If Ln[smg](q) • n e < 0, then d/dt (D e (q + t Ln[smQ](<7))) 4=[) = n e ■ Ln[smQ](q) < 0, 
which implies that D e (g + tLn[smg](g)) < D e (q) = sm Q(<z) for t > small enough. 
On the other hand, invoking Theorem 12.61 for the function — sitlq, we have that 
sm.Q(q) + t\\ Ln[smQ](g)|| 2 /2 < sm Q (q + t Ln[smQ](g)) < B e (q + t Ln[smQ](g)). Gath- 
ering both facts, we conclude t\\ Ln[sniQ] (q) || 2 /2 < 0, which is a contradiction. 

Let us now prove (iv). By definition, D e (<?) < IR(Q). This inequality is strict if 
q g" IC(Q). Let x <E IC(Q). If we take a point O in the edge e, then the function D e 
can be expressed as D e (p) = (p — O) ■ n e . Then, we have 

D e (x) = (x — O) -n e = (x — q) ■ n e + (q — O) -n e = {x — q) ■ n e + B e (q) . 

Since D e (a;) > smg(x) = IR(Q), we conclude that (x — q) ■ n e > IR(Q) — D e (g) > 0, 
and that the second inequality is strict if q g" IC(<5)- □ 

4. Analysis of the multi-center functions. Here we study the locational op- 
timization functions Hdc an d Ti-sp for the disk-covering and sphere-packing problems. 
We characterize their smoothness properties, generalized gradients, and critical points 
for arbitrary numbers of generators. 

4.1. Smoothness and generalized gradients. We start by providing some 
alternative expressions and useful quantities. We write 

H DC (P)= max G t (P) , H S p(P) = min F t (P) , 
<e{i,. ..,«.} ie{i,...,n} 

where 

Gi(P)= max Fi(P)= min \\ q - Pi \\ . 

q€Vi(P) q£mt(Vi(P)) 

Note that Gi(P) = lg Vi (p)(Pi) an d Fi{P) = sm 7 .( P )(pi), where, for i e {1, . . . ,rt}, 

lgvr 4 :Vi->R, sm^-.V^R. 

Proposition ^ . 3l providcs an explicit expression for the generalized gradients of lg v . and 
smy; when the Voronoi cell Vi is held fixed. Despite the slight abuse of notation, it is 
convenient to let d\g Vi ^ P ^(pi) denote dlg v (jpi)\v=v i ( L P)i an( i let d sm Vi(P)(Pi) denote 
dsmv(Pi)\v=Vi(p)- 

In contrast to this analysis at fixed Voronoi partition, the properties of the func- 
tions Gi and F{ are strongly affected by the dependence on the Voronoi partition 
V(P). Wc endeavor to characterize these properties in order to study Tiuc an( l Ti-sp. 

PROPOSITION 4.1. The functions Gi,F t : Q n — > M are locally Lipschitz and 
regular. As a consequence, the locational optimization functions TLdc^SP '■ Q" "~ * K 
are locally Lipschitz and regular. 

Proof, (a) Gi is locally Lipschitz and regular. The definition of the function Gi 
admits the following alternative expression 

G i {P)= max \\pi-v\\ . (4.1) 

v£Vc(Vi) 
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Let Pq be nondegenerate at the ith generator. Then, there exists a neighborhood 
U of Pq where the set Af(i) does not change. Let . . . , Dj^J, {u>i, . . . , wm 2 }-, 
{z%, . . . ,Zm 3 } be the vertexes of Vi of types (a), (b) and (c) respectively. Then, Gi 
can be locally written as 



Gi(P)=max^ max max max \\zt-p, 

{££{!,. ..,Mi} te{l,...,M2} "te{l,...,M 3 } 

for all P G U. Therefore, Gi restricted to U coincides with the function GuU) '■ Q r 
M. defined by 

m {P) = max| f6 max Mi} \\v t - Pi \\, ^rnax^ \\w t - M, u *»*, a} \\» ~ 

(4.2) 

The function Gnu) is the maximum of a fixed finite set of locally Lipschitz and regular 
functions, and consequently, locally Lipschitz and regular by Proposition 12.41 We 
conclude that Gi is both locally Lipschitz and regular at Po- 

Let Pq be degenerate at the zth generator. Then, in any neighborhood U of Pq 
there are different sets of neighbors of the ith generator. Indeed, because the number 
of generators, edges of the boundary Q and vortexes of Q is finite, there is only a finite 
number of different sets of neighbors of the ith generator over U, say Af 1 (i), . . . ,J\f L (i). 
This implies that Gi admits the following alternative expression over U 

Gi(P)=min{g MHi) (P),...,g^ L(i) (P)} . (4.3) 

Again resorting to Proposition EH we conclude that Gi is both locally Lipschitz and 
regular at Pq. 

(b) Fi is locally Lipschitz and regular. From the definition of Fi, it is clear that 
its value at a configuration P is attained at the boundary of the Voronoi region Vi. 
Therefore, one only minimizes among the edges associated with the Voronoi neighbors 
Af(i) and the edges of Q with non-empty intersection with V{. Moreover, one can also 
see that the minimum must be attained at a point of the form proj e (pi), for some 
edge e of Vi- Now, consider the function J^^fu) ■ Q n — ► R defined by 



(P) = min I min \\ Pl - - - Pl || , min D e (pi) \ 



(4.4) 



We shall prove that ^-jv(i) coincides with Fi. If k £ M(i), then (pi +pf-)/2 £ Vi. 
Since Q \ Vi is open, there exists a neighborhood of (pi + pk)/2 such that U C Q/V.. 
Therefore, 

_ £l+£i|| > ^ \\ p . _ g|| > i _ = Fi{P) . 

I qeU <}0mt(Vi) 

If an edge e of Q does not intersect Vi, then proj e (pi) ^ Vi. Using again the fact that 
Q\Vi is open, there exists a neighborhood U of proj e (pi) in M 2 such that UPiQ C Q\Vi. 
Then, 

lift - Proj e (p 4 )|| > min \\ Pi - q\\ > F^P) . 
qeunQ 

As a consequence of the previous inequalities, Fi equals Tj^^y Being the minimum 
of a fixed finite number of locally Lipschitz and regular functions on Q n , Fi is also 
locally Lipschitz and regular by Proposition ^. 41 □ 



16 



Jorge Cortes and Francesco Bullo 



Next, one can actually prove the following stronger result. 

Proposition 4.2. The locational optimization functions Hd&Usp '■ Q n "~ * ^ 
are globally Lipschitz, with Lipschitz constant equal to 1. 

Proof, (a) TLdc * s globally Lipschitz. Let P, P' be two configurations of the n 
generators. Without loss of generality, assume that TLdc{P) < Wdc(-P')- Let i, j 
and g , q' G Q be such that W D c(P) = = ||go ~Pi|| and Wdc(P') = Gj(P') = 

Iko — Pj'll- Now, consider the set B-2(q' ,Gi(P)). Then there exists a /c such that 
Pk G P>2(q' a ,Gi(P)) (otherwise, ||g — > Gj(P), which contradicts the definition of 
the function 7Ydc)- On the other hand, we necessarily have that p' k B2(q' ,Gj(P')), 
since otherwise \\q' Q — p' k \\ < \\q' — p'j\\, which implies that q' £" Vj, contradiction. 
Finally, we apply the triangle inequality to obtain ||qq — p' k \\ < \\q' Q — pf.\\ + \\pk — pl||- 
Gathering the previous facts, we have 

IWdc(P') - Wdc(P)| = - Gi(P) 

< ll?o -p'd U -PkW < \\Pk -p' k \\ < \\p- p'W ■ 

(b) TLsp is globally Lipschitz. Let P, P' be two configurations of the n generators. 
Without loss of generality, assume that Hsp(P) < Ttsp(P')- Let i, j and go, Qq G Q 
be such that H SP (P) = F 4 (P) = ||g and H SP (P') - F,-(P') = ||g We 

treat separately the following two cases: (i) qo does not belong to the boundary of Q, 
and (ii) go belongs to the boundary of Q. In case (i), it necessarily exists k £ Af(i) 
such that ||go = \\qo-Pk\\. If \\qo -p-|| > F,(P'), then 

IWsp(P') - W SP (P)| = Fj(P') - Fi(P) < ho-p'iW - \\qo-Pi\\ 

<\\Pi-p^\\<\\P-P'\\. (4.5) 

If, on the contrary, ||g — Pill < Fj(P'), then g £ mt(t^'). Therefore, ||g — p' k \\ > 
Pfc(P') > Pj(P')- Now, we perform the same computation as in l|4.5[) to conclude 
|«sp(P')-«sp(P)| < ||P-P'||. 

In case (ii), we prove that ||go — Pill > Fj(P'). Suppose this is not true, i.e., 
1 1 So — Pill < Fj{P')- Let m = g + e(go — Pi)j with sufficiently small e > such that 
||m-p-|| < Fj(P'). Clearly m £ Q. On the other hand, by definition B 2 {p' i , F t {P')) C 
Vj'. Now, we have, 

B 2 (p' l ,F 3 (P')) c B 2 (p' l ,F l (P')) CV- CQ. 

But, since ||m — p' || < Fj(P'), then to € B 2 {p'i, Fj(P')) C Q, which is a contradiction. 
Therefore, ||go — Pi\\ > Fj(P'), and now the same argument as in (|4.5(l guarantees 

that \n SP (P')-n SP (p)\ < ||P-P'||. 

We now introduce some quantities that arc useful in characterizing the generalized 
gradient of the functions Gi. Given a vertex of type (b), v = v(e,i,j), determined 
by the edge e and two generators pi and pj, we consider the scalar function \{e,i,j) 
defined by 

P r °j e (Pi ~ v(e,i,j)) = X(e,i,j) proj e (pj - pi) (4.6) 

where P e is the orthogonal projection onto the edge e; see Fig. 14.11 One can see that 
A(e, i, j) + A(e, j, i) = 1. If e is a segment in the line ax + by + c = 0, (Axij , Ay^) = 
Pj —p^ (x m ,y m ) = (pi +pj)/2, then one can show 

. _ 1 _ (aAxjj + bAy lj )(ax m + fa/ m + c) 
le ' Z '- ?j_ 2 (aA yij - bAxij) 2 
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Fig. 4.1. To illustrate eq. 14.61 we draw the vectors proj e (pj — v(e,i,j)) and proj e (pj — Pi) 
for various locations ofpi, pj, and e. The left, center and right figures correspond to \(e,i,j) > 0, 
\{e,i, j) = 0, X(e,i,j) < 0, respectively. 



Given a vertex of type (a), v = v(i,j,k), determined by the three generators Pi, pj, 
and p k , we consider the scalar function fi(i,j, k) defined by 

P r °je 3fc (Pt ~ v (hj> k )) = KhJ> k ) P r °je Jfc (Pi - Pi) 

where ejk is the bisector of pj and pk and where pi = pj if pj belongs to the half-plane 
defined by ejk containing pi, and pi = pk otherwise. One can see that fi(i,j,k) = 
n(i, k,j) and that /J,(i,j, k) + k, i) + n(k, i,j) = 1. From the expression for A, one 
can obtain 

fiU j k) = -+ (^ Xi i Ax i k + AyijAyjk)(Axi k Ax jk + Ay ik Ay jk ) 
2 2(x k Ayij - XjAy ik + XiAy ]k ) 2 

Note that, in general. A and \x are not positive functions. Now we are ready to describe 
in detail the structure of the generalized gradient of the functions G,;, Fj. 

Proposition 4.3. The generalized gradient of Gi : Q n — >• R at P e Q n is 

dG l {P) = co {d v Gi(P) G (R 2 )™ I « G Ve(V t (P)) such thatG l (P) = \\ Pi -v\\} 

where we consider separately the following cases. If V = v(i,j,k) is a nondegenerate 
vertex of type (a), then 

d v (i,j,k)Gi(P) = dv(k,i,j)Gk(P) = dv( 3 .k,i)Gj(P) = 
(0, . . .,fi(i,j,k)vTs{pi -«),.. .,fi(j,k,i)vTs{pj -v),.. . ,fi(k,i, j) vrs(p fe -«),..., 0) 

S v ' * * ' S v ' 

ith place jth place kth place 

where, without loss of generality, we let i < j < k. If v = v(e,i,j) is a nondegenerate 
vertex of type (b), then 

9 v ( e ,i,j)Gi{P) — dv{e,j,i)Gj(P) 

= (0,...,\(e,i,j)vrs(pi - v), . . . , A(e, j, i) vrs(pj -v),...,0) 

S ^ i ^ 

ith place jth place 

where, without loss of generality, we let i < j. If v = v(e,f,i) is a nondegenerate 
vertex of type (c), then 

d v(eJA Gi{P) = (0,...,0 > vrs(pi-«) ) 0,...,0). 

ith place 

Finally, if the vertex v is degenerate, i.e., if v is determined by d > 3 elements 
(generators or edges), then there are (f 1 ^ 1 ) P a i rs °f elements which determine the 
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vertex v together with the generator p^. In this case, d v Gi(P) is the convex hull of 
d v {a,f3.~t)Gi{P) for all (''j 1 ) such triplets (a,/?, 7). 

Note that, at all nondegenerate configurations P, the quantity d v Gi{P) is the 
generalized gradient of the function (pi,...,p n ) 1— > \\pi — v(i,j, k)\\; however, this 
interpretation cannot be given when P is degenerate. 

Proof. We present the proof for the expression for dGi(P). Let us consider first 
the case when P is nondegenerate configuration for the zth generator. According to 
the proof of Proposition ^. II Gi coincides with the function Gu(i) over a neighborhood 
U oiP. Hence, dGi(P) = dGj\f(i)(P) which, according to eq. (|4.2|) and Proposition l2T4l 
takes the form 

■£s\\v-pi\\ I v g Ve{Vt(P)) such that ||«-p<|| = G^P] 
oP 

If v = v(i,j, k) is a nondegenerate vertex of type (a), then 
d ( dv \ 

g-\\Pi ~ v (h3> k )\\ = vrs fe ~ v ) [ J 2 - g^-J = V(ij,k)vrs(pi - v) , 

d ( dv \ 

Q^WPi _ v ( l J> k )\\ = -vrsfe - v) I — I = fj,(j,k,i)vrs{pj - v) , 

d 

-z—\\Pi-v{i,j,k)\\=Q, £^i,j,k. 
ope 

If v = v(e, is a nondegenerate vertex of type (b), then 
d ( dv \ 

g—\\Pi - v(e,i,j)\\ = vrs(pi ~ v) \ I 2 - g-\ = X{e,i,j)vTs{Pi - v) , 

d ( dv \ 

\\Pi ~ v(e,i,j)\\ = - vrs(p, -?;) — = \(e,j,i)vTs{pj - v) , 



dp, V^ft 

— \\pi -v{e,i,j)\\ = 0, 
apt 

If v = v(e, /, i) is a nondegenerate vertex of type (c), then 

— - t)(e, /, i) || = vrs(p 4 - v) , 

^j||j>i-«(e I / J i)|| = > 

If P is a degenerate configuration at the ith generator, then this function can 
be expressed as in cq. (|4.3|l in a sufficiently small neighborhood U of P. According 
to Proposition 12.41 the generalized gradient of Gi is given by the convex hull of 
the generalized gradients of each of the functions Gm 1 m , . . . , Gj^l . The claim now 
follows by reproducing the previous discussion for the generalized gradients of each of 
the functions Gj^i(i), £ G {1, . . . , L}. □ 

The expression for dFi(P) can be deduced in an analogous (and simpler) way, 
since according to the proof of Proposition ^. II it is not necessary to establish any dis- 
tinction between the degenerate and the nondegenerate configurations. Accordingly, 
we state the following result without proof. 

PROPOSITION 4.4. The generalized gradient of P, : Q n — > R at P e Q n is 

dFi(P) = co{d e Fi(P) e (K 2 )™ I e e Ed(^(P)) such that F^P) = B e ( Pl )} 
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where, if e = e{i,j) is an edge of type (a), then 

d e (z,j)Fi(P) = 9 e ( jti )Fj(P) = i(0, . . . , n e (i tj ) -n e( i d ),. . . ,0), 

ith place jth place 

and if e = e(i) is an edge of type (b), then 

d e (i)Fi(P) = (0, . .., n e(i) , . . . ,0). 

ith place 



Next, we give conditions under which the functions A and /i take positive values. 
Lemma 4.5. Let P G Q" and let v G Ve DC (V(P)). Then, 

(i) if v belongs to an edge e ofQ, then there exist generators pi and pj such that 
X(e,i,j) and\(e,j,i) are positive, and 

(ii) if v belongs to int(Q), then there exist generators pi, pj and pk such that 
/j,(i,j,k), n(j,k,i) and fi(k,i,j) are positive. 

Proof. Consider first the case when v is nondegenerate. If v is in the edge e of 
Q (i.e. v is of type (b)), let pi and pj the two generators determining it. From the 
definition of A. one sees that the values A(e,i, j) = and A(e,j, i) = correspond to, 
respectively, pj and pi lying on the orthogonal line to e passing through v(e,i,j). If 
A(e,i, j) < 0, then there exists w G e n Vj such that \\pj — w\\ > \\pj — v\\ = Hnc(P), 
which is a contradiction. Therefore A(e,i,j) > 0. The same argument guarantees 
A(e,j,z) > 0. If v is of type (a), and Pi, pj and pk are the elements determining it, 
a similar argument leads to the conclusion that fc), k,i) and fi(k,i,j) are 

positive. 

Consider the case when v is degenerate. Let . . . , i m } be such that v € VL, 
j G {1, . . . , m}. Assume v is an edge e of Q. Let I denote the orthogonal line to the 
edge e passing through v. We claim that there must exist generators in . . . ,i m } 
on both sides of I. Assume this is not the case, i.e., {pi t , ■ ■ ■ ,Pi m } are contained in one 
of the closed half-planes defined by I, say i_. Take w G Z+ P\ e arbitrarily close to v. 
Since {pi t , ■ ■ ■ ,Pi m } C 1-, we have \\pi j —w\\ > \\pi j —v\\. On the other hand, since no 
generator outside the set {pi tl ■ ■ ■ ,Pi m } is involved in the definition of v, there must 
exist j* such that w G V t ., . Therefore, G tj , (P) > \\p ijt —w\\ > \\p ijt -v\\ = Hbc(P), 
which is a contradiction. Assume now that v G int(Q). Our claim is that, for any 
line / passing through v, there must exist generators on both sides of I. If this is not 
the case, i.e., {pi t , ■ ■ ■ ,Pi m } C I—, then take w G B^{v, e) fl Z+ fl o, where o denotes 
the orthogonal line to I passing through v. As before, w G Vi., for some j* and 
\\pi-, — w|| > \\pij, — v\\, which yields a contradiction. □ 

This completes our analysis of the generalized gradients of Gi and Fi and, with 
these results, we return to studying the generalized gradients of Hdc and Wsp- An 
immediate consequence of Propositions 12.41 and 14. II is that 

dHvc(P)=co{dGi(P)\i£l(P)} , 

dH SP (P) - co {dFi(P) | i G I(P)} . (4.7) 

Furthermore we can provide the following more detailed characterization. 

PROPOSITION 4.6. Let P G Q n . For each i G {1, . . . , n}, the image by m of the 
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generalized gradients ofTloc an d Hsp °^ P * s given by 

( n (8Gi(P)) if i £ I(P) , Ve DC (V(P)) C Ve(Vi{P)) 

MdH DC (P)) = < co{7Ti(dGi(P)),0} if i e I(P) , Vc DC (V(P)) <jL Ve(Vi(P)) 

[O ifi^I(P) 

f TTi(dFi(P)) if i £ I(P) , Ed SP (V(P)) c m(Vi(P)) 

^(dU SP {P)) = I coin(dFi(P)),0} if i e I(P) , Ed SP (V(P)) <jL Ed(V t (P)) 

[0 ifi^I(P) 

Proof From cq. |HJ, if i £ I(P), then 7ri(5W D c(P)) = 0, 7r 4 (cW SP (P)) = 0. If 
i £ P\P), then using Proposition ^. 31 we deduce that the generators^ such that dGj 
has a nonzero entry in the ith place (and hence contributes to the projection by iti 
of dHr>c) must share a vertex with the ith generator. Analogously, if i £ I(P), then 
using Proposition ^. 41 we deduce that the generators pj such that dFj has a nonzero 
entry in the ith place (and hence contributes to the projection by ~Ki of dTtsp) must 
satisfy j £ J\f(i). For the disk-covering function, if v is a common vertex of Vi and 
Vj, determined by i, j and a third element a, then d v ( a j^Gj = d v ( a ^^Gi, and 
the expression for TTi(dTtr>c(P)) then follows. The argument for the expression of 
TTi(dHsp(P)) is analogous. □ 

4.2. Critical points. Having characterized the generalized gradients of Hbc 
and Hsp, we now turn to studying their critical points. 

Theorem 4.7 (Minima of Hue)- Let P £ Q n be a nondegenerate configuration 
and £ mt (dH.Dc{P))- Then, P is a strict local minimum ofTLoc^ a ^ generators 
are active and P is a circumcenter Voronoi configuration. 

Proof. Since P is nondegenerate, note from Proposition 14.31 that d v Gi is a sin- 
gleton for each v £ Ve(Vi(P)), i £ {1, . . . ,n}. Let w £ (K 2 )". We claim that moving 
the configuration of the generators from P in the direction w can only increase the 
cost. The hypothesis £ mt(dH.Bc{P)) implies by Lemma [2. II that there exists i 
and v £ Ve(Vi(P)) n Vc DC (V(P)) such that w ■ d v G l (P) > 0. Since P is nondegen- 
erate, v will still belong to Vi(P + ew) for sufficiently small e > 0, and consequently 
Ht>c{P + ew) > Gi(P + ew) > G,(P) = Hnc(P)- Therefore P is a strict local 
minimum. 

Since Tti is an open map, the set 7Ti(int(cWDc(P))) is open for each i £ {1, . . . , n}. 
Therefore, 7T,(int(97iDc(P))) 0> an d hence all generators are active, i.e. I(P) = 
{1, . . . , ??}. Let us see that all generators must also be centered. Assume P is non- 
degenerate and consider the ith generator. Take w £ M 2 and let w £ (R 2 )™ be the 
vector with has w in the ith place and otherwise. By Lemma 12. II there exist j and 
v £ Ve(Vj(P))nVe D c(V(P)) such that w -d v Gj > 0. Since W-d v Gj = w-Ki{d v G-) > 0, 
then necessarily TTi(d v Gj) ^ 0, and therefore v £ Vi(P) and ni(d v Gj) = TVi(d v Gi). The 
vertex v is determined by pi, pj and a third clement, say a. Depending on whether 
a corresponds to an edge or to another generator, we have that iri(d v Gi) is equal to 
\(a, i,j) vrs(pi —v) or fi(a, i,j) vrs(jpi — v). In any case, from Lemma T4.5I we deduce 
that X(a,i,j) (respectively fi(a,i,j)) belongs to (0,1). Therefore w ■ TTi(d v Gi) > 
implies w ■ vrs(pi — v) > 0. Consequently, £ mt(dlg v ./ P \(pi)). By Proposition 13.31 
this implies that pi = CC(V^). Hence, P is a circumcenter Voronoi configuration. □ 

Theorem 4.8 (Maxima of H S p). Let P £ Q n and £ mt(dH S p(P)). Then, P 
is a strict local maximum ofTLsp, a ^ generators are active and P is a generic incenter 
Voronoi configuration. 
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Fig. 4.2. Local extrema of the disk- covering and the sphere-packing functions in a convex 
polygonal environment. The configuration on the left corresponds to a local minimum of Hoc with 
£ dTi.r>c{P) an d m t{dHnc{P)) = The configuration on the right corresponds to a local maxi- 
mum ofHgp with £ dHsp(P) and mt(&Hsp(P)) = 0- In both configurations, the 4th generator is 
inactive and non- centered. 



Proof. The proof of this result is analogous to the proof of Theorem 14.71 Note 
that G mt(9smy.(p)(pj)) implies, by Proposition 13.31 that lC(Vi(P)) = {pi}, and 
hence P is a generic incenter Voronoi configuration. □ 

Remark 4.9. Theorems 14.71 and 14.81 precisely provide the interpretation of the 
multicenter problems that we gave in Section f2.2l since all generators are active, they 
share the same radius. If one drops the hypothesis that belongs to the generalized 
gradient of the locational optimization function, then one can think of simple examples 
where P is a local minimum of Hdc (respectively local maximum of Hsf), and there 
are generators which are inactive and non-centered, see Fig. 14.21 

5. Dynamical systems for the multi-center problems. In this section, we 
describe three algorithms that (locally) extremize the multi-center functions for the 
disk-covering and the sphere-packing problems. We first examine the gradient flow 
descent associated with the locational optimization functions Wdc an d Hsp. This flow 
is guaranteed to find a local critical point, but it has the drawback of being centralized, 
as we describe later. Then, we propose two decentralized flows for each problem. One 
roughly consists of a distributed implementation of the gradient descent. As we show, 
it is very much in the spirit of behavior-based robotics. The other one follows the 
logical strategy given the result in Theorems l4.7l and l4.8l each generator moves toward 
the circumcenter (alternatively, incenter set) of its own Voronoi polygon. We call them 
Lloyd flows, since they resemble the original Lloyd algorithm for vector quantization 
problems, where each quantizer moves toward the centroid or center of mass of its 
own Voronoi region, see ^] 12] < We present continuous-time versions of the 
algorithms and discuss their convergence properties. In our setting, the generators' 
location obeys a first order dynamical behavior described by 

pi = Ui(pi, . . . ,p n ) , ie{l,...,n}. (5.1) 

The dynamical system <|5.1ll is said to be (strongly) centralized if there exists at 
least an i G {1, . . . , n} such that Ui(j>i, . . . ,p n ) cannot be written as a function of 
the form Uiijp^p^, . . . ,Pi m ), with m < n — 1. The dynamical system (|5.1|) is said 
to be Voronoi- distributed if each . . . ,p n ) can be written as a function of the 

form Ui(pi,Pi 1 , . . . ,pi m ), with i^ G N{P,i), k G {1, . . . , m}. Finally, the dynamical 
system (|5.1|) is said to be nearest-neighbor- distributed if each . . . ,p n ) can be 

written as a function of the form Ui^pi^p^, . . . ,Pi m ), with \\pi ~Pi k \ \ < \\pi —Pj\\ for all 
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j G {1, . . . , n}, and k 6 {1, . . . , m}. A nearest-neighbor-distributcd dynamical system 
is also Voronoi-distributed. 

It is well known that there are at most in — 6 neighborhood relationships in 
a planar Voronoi diagram |23l see Section 2.3]. Therefore, the number of Voronoi 
neighbors of each site is on average less than or equal to 6. (Recall that sites are 
Voronoi-ncighbors if they share an edge, not just a vertex.) We refer to for more 
details on the distributed character of Voronoi neighborhood relationships. 

Note that the set of indexes {ii, ■ ■ ■ ,i m } f° r an specific generator pi of a Voronoi- 
distributed or a nearest-neighbor-distributed dynamical system is not the same for all 
possible configurations P. In other words, the identity of both the Voronoi neighbors 
and the nearest neighbors might change along the evolution, i.e., the topology of the 
dynamical system is dynamic. 

5.1. Nonsmooth gradient dynamical systems. Consider the (signed) gen- 
eralized gradient descent flow 1|2.7[) for the locational optimization functions 7Ydc and 

P=-Ln(cW DC )(P), P = Ln(0W SP )(P) . 

Alternatively, we may write for each i £ {1, . . . , n}, 

Pi = -7Ti(Ln(#W D c)(pi, ■ ■ ■ ,Pn)) , (5-2) 
pi = n i (hn(dHsp)(pi 1 ■ ■ ■ ,p n )) ■ (5.3) 

As noted in Sect ion ITU these vector fields are discontinuous, and therefore their solu- 
tion must be understood in the Filippov sense. Eq. I|4.7|l and Propositions 14. 31 and !4. 41 
provide an expression of the generalized gradients at P, dTlp,c{P) an d dHsp(P)- 
One needs to first compute the generalized gradient, then compute the least-norm 
element, and finally project to each of the n components; therefore the expressions in 
Proposition ^. 61 arc not helpful. Note that the least-norm element of convex sets can 
be computed efficiently, see [H], however closed- form expressions are not available in 
general. 

One can see that the compact set Q n is strongly invariant for both vector fields 
— Lii(cWdc) and Ln(cWsp)- Regarding — Ln(cWoc), this is a consequence of Propo- 
sition l4~3l and of Lemma [4. 51 Regarding Ln(dTisp), this is a consequence of Proposi- 
tion E31 

Proposition 5.1. For the dynamical system (|5.2() (respectively (|5.3|) ). the gen- 
erators' location P ~ (pi, . . . ,p n ) converges asymptotically to the set of critical points 
ofTi-DC (respectively, ofH S p)- 

Proof. From Propositions 14. II and 14.21 we know that Hbc an d T~Csp arc globally 
Lipschitz and regular over Q n . The result follows from Proposition 12 . 1 fl considering 
the dynamical system restricted to the strongly invariant and compact domain Q n . 
□ 

Remark 5.2. The gradient dynamical systems enjoy convergence guarantees, but 
their implementation is centralized because of two reasons. First, all functions Gi{P) 
(respectively Fi(P)) need to be compared in order to determine which generator 
is active. Second, the least-norm element of the generalized gradients depends on 
the relative position of the active generators with respect to each other and to the 
environment. 

Remark 5.3. As illustrated in Fig. 15. li the evolution of the gradient dynamical 
systems may not leave fixed even the generators that are centers (circumcenter or 
inccnters). 
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Fig. 5.1. Illustration of the gradient descent. In the left figure, the jth generator is in the 
circumcenter of its own Voronoi region, but the control law 15.21 drives it toward the vertex v. 
In the right figure, the jth generator is in the incenter of its own Voronoi region, but the control 
law 15.31 drives it away from the edge e. 

5.2. Nonsmooth dynamical systems based on distributed gradients. In 

this section, we propose a distributed implementation of the previous gradient dy- 
namical systems and explore its relation with behavior-based rules in multiple- vehicle 
coordination. Consider the following modifications of the gradient dynamical sys- 
tems (jjT2l - (jOt . 



is determined only by the position of pi and of its Voronoi neighbors Af(P, i). On the 
other hand, the system H5.5|) is nearest-neighbor-distributed, since Ln(9smy j (p))(P) 
is determined only by the position of pi and its nearest neighbors. 

For future reference, let Ln(<91g v )(P) = (Ln(91gy i(P) )(P), . . . , Ln(51g Vn(P) )(P)), 
Ln(9smy)(P) = (Ln(9sm V ' 1 (_p))(P), . . . , Ln(9smy„(p))(P)), and write 



As for the previous dynamical systems, note that these vector fields are discontinuous, 
and therefore their solutions must be understood in the Filippov sense. One can 
see that the compact set Q n is strongly invariant for both vector fields — Ln(dlgy) 
and Ln(9smy). This fact is a consequence of the expressions for the generalized 
gradients of lg and sm in Proposition 13.31 Note that in the 1-center case, (|5.2|l 
(respectively ()5.3fl ) coincides with l|5.4|l (respectively with l|5.5|l l. 

Proposition 5.4. Let P e Q". Then the solutions of the dynamical sys- 
tems (|5.4|) and (I5.5|l starting at P are unique. 

Proof, (a) Uniqueness of solution for l|5.4|l . Let D\ g be the set of P £ Q n such that 
P is nondegenerate and lgy (p)(.Pi) is attained at a single vertex for all i. Note that 
Q n \D\ g has measure zero, and that the vector field — Ln(<91g v ) is differentiable (and 
hence locally Lipschitz) when restricted to any connected component of D\ g . Let P, 
P' belong to different connected components of D\ g , and let [|P — P'|| < e. Consider 
all the indexes i at which the values of lgy(p)(Pi) and \g,Vi(P')(Pi) are attained at 
different vertexes. For these indexes, 




P = -Ln(<91g v )(P), P = Ln(<9sm v )(P). 



Ln(91gy. ( p))Oi) + Ln(dlgy. ( p,))(p-) 



VTS(V - Pi) - VTS(W' - p[) , 
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for certain vertexes v and w' . Note that for e small enough, the vertex w' in the 
Voronoi configuration P' corresponds to a vertex w in the Voronoi configuration P. By 
construction, pi and p[ belong to an 0(e) neighborhood of the bisector b vw determined 
by v and w, and n vw ■ (jpi — p[) < 0. In addition, the component of vrs(u — pi) — 
vrs(u/ — p'j) along b vw is O(e) whereas n vw ■ vrs(v — pi) > and n vw ■ vrs(«/ — p[) — 
n vw ■ vrs(w — pi) + O(e), with n vw ■ vrs(«; — pi) < 0. Then, 

vts(v — pi) — vrs(w — p[) 

= proj„ um (vrs(u — pi)— vrs(V - p[) + proj biu (vrs(i> - p,) - vrs(w' - p[) 

= proj nijm (vrs(i> - pi) - vrs(V - p\) + 0(e) , 

and, in turn, for sufficiently small e 

(Pi - Pi) ■ ( v rs(v - Pi) - vrsfV - p-)) 

= (n vw ■ (pi - p'i)){n vw ■ (vrs(w - p f ) - vrs(u/ - p.-))) + 0(e 2 ) < . 

The result now follows from Theorem 1 at page 106 in |15| . 

(b) Uniqueness of solution for (|5.5|) . Let D sm be the set of P £ Q n such that 
smy.(p)(j)i) is attained at a single edge for all i. Note that Q n \ D sm has measure 
zero, and that the vector field Ln(<9smy) is differentiable(and hence locally Lipschitz) 
when restricted to any connected component of D sm . Let P. P' belong to different 
connected components of £> sm , and let \\P — P'\\ < e. Consider all the indexes i at 
which the values of smy.f P \(pi) and smy.(p/)(p^) are attained at different edges. As- 
sume these edges are of type (a) (the type (b) case can be treated analogously). For 
these indexes, 

Ln(dsm Vi (p))(pi) - Ln(dsm Vi{P , ) )(p / l ) = vrs(p l - pj) - vrs(p- -p' k ) , 

for some uniquely determined pj and p' k , with j ^ k. By construction, pi and p[ 
belong to an O(e) neighborhood of the bisector bjk determined by pj and pk, and 
n kj • {Pi — Pi) < 0. In addition, the component of vrs(p^ — pj) — vrs(pj — p'k) along bjk 
is O(e) whereas n kj ■ vrs(p, —pj) > and n k j ■ vrs(p- -p' k ) = n kj ■ vrs(p; -p k ) +0(e), 
with n k j ■ vrs(pi — p k ) < 0. Then, 

vrs(p, - Pj)- vrs(p- - p' k ) 

= Pi'°j„ fcJ ( vrs (Pi - Pj) ~ vrs(p- - p' fc ) + proj bjfc (vrs( Pl - ft) - vrs(p- - p' fe ) 

= Pr°jr ifcJ (vrs(pi -Pj) - vrs(p* - p k ) + O(e) , 

and, in turn, for sufficiently small e 

(Pi - Pi) ■ (vrs( Pi - pj) - vrs(p- - p' fc )) 

= (ny ■ (pi - p'i)){n kl ■ (vrs(pi - p-,) - vrs(p, - p fe ))) + 0(e 2 ) < . 

The result now follows from Theorem 1 at page 106 in ^3- D 

Remark 5.5 (Relation with behavior-based robotics: move toward the furthest- 
away vertex) . The distributed gradient control law in the disk-covering setting (|5.4J) 
has an interesting interpretation in the context of behavior-based robotics. Consider 
the ith generator. If the maximum of lg y-^p) is attained at a single vertex v of its 
Voronoi cell Vi, then lgy.^p) is diffcrentiable at that configuration, and its derivative 
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corresponds to vrs(pi— v). Therefore, the control law l|5.4|l corresponds to the behavior 
"move toward the furthest vertex in own Voronoi cell." If there arc two or more 
vertexes of Vi where the value lgv; (p)(Pi) is attained, then (|5.4|) provides an average 
behavior by computing the least-norm element in the convex hull of all vrs(j»i — v) 
such that \\pi - v\\ = lg v . (P) (pj). 

Remark 5.6 (Relation with behavior-based robotics: move away from the nearest 
neighbor) . The distributed gradient control law in the sphere-packing setting i|5.5[) has 
also an interesting interpretation. For the zth generator, if the minimum of smy.(p) 
is attained at a single edge e, then smy.(p) is diffcrcntiable at that configuration, 
and its derivative is n e . The control law l|5.5Jl corresponds to the behavior "move 
away from the nearest neighbor" (where a neighbor can also be the boundary of the 
environment). If there are two or more edges where the value smy.(p)(pi) is attained, 
then H5.5f) provides an average behavior in an analogous manner as before. 

Proposition 5.7. For the dynamical system l|5.4[) . the generators' location P = 
(pi, . . . ,p n ) converges asymptotically to the largest weakly invariant set contained in 
the closure of Av C {Q)^= {P eQ n \ie I(P) => Pl = CC(V<)}- 

Proof. Let a G £-Ln(dig v )WDc(P)- By definition, a = - Ln(<9 lg v )(P) • C, for 
all C G dHuc(P)- Let v G Ve D c(V(P)). From Lemmas PI and H751 we know that, 
independently of the dcgcncratc/nondcgcnerate character of the Voronoi partition at 
v, there always exist either an edge e of Q and generators pi and pj , or generators pi, pj 
andpfc, such that X(e,i,j), X(e,j,i) > (respectively fx(i,j,k), p(j,k,i), p(k,i,j) > 
0). If v is a vertex of type (b), then 

a=-Ln(d\g v )(P)-d v G l (5.6) 
= -Ln(aig Vi(P ))(P) ■ X(e,i,j)vrs(pi - v) - Ln(dlg Vj{P) )(P) ■ X(e,j,i) vrsfe -v). 

From Lemma I3.8f i1 wc conclude that a < 0, and the inequality is strict if either 
Pi ^ CC(Vi) or pj ^ CC(Vj). The same conclusion can be derived if v is a vertex 

of type (a). Therefore, max/L L n(ai gv )^Dc(P) < or jC_ hn{aigv) H D c{P) = 0- 
Now, resorting to the LaSalle principle ( Theorem 12. 9|l . wc deduce that the solution 
P : [0, +oo) — > Q n starting from Pq converges to the largest weakly invariant set 
contained in Z_ Ln(aig v ),w D c 

Let us see that Z_ Ln(aig v ),w DC n Q n ^ s ec l ua l to Ar>c(Q)- Take a configuration 
P G A DC (Q). Then, Ln(dlg Vi(P) )(P) = if * G J(P), and tt^C) = if * ^ I(P), for 
any ( G dHr>c{P) (cf. Proposition ^. 6|1 . Consequently, = — Ln(51g v )(P) • C, for all 
C G dHvc(P), andsoO G Ln(aigv) H D c(P)- Therefore, A DC (Q) C ^- L n(aig v ),w DC - 
Now, consider P G Z_ Ln(dlgv y HDC . Then, G £-Ln(ai gv )^Dc(P), that is, = 
— Ln(<91g v )(P) • C, for all C G cWdc(P)- If P is nondegenerate, we deduce from 
eq. (|5.6|l and Lcmma lTHl that all the active generators are centered, i.e., P G ^dc(Q)- 
If P is degenerate, consider a degenerate vertex v where the value of Hbc{P) is 
attained. For simplicity, we deal with the case where v is contained in an edge e 
of Q (the case v G int(Q) is treated analogously). From Lemma [4.51 we know that 
there exist generators pi, pj determining v on opposite sides of /, the orthogonal line 
to the edge e passing through v. From eq. (|5.6() and Lemma 13.81 we deduce that 
both pi and pj are centered. Now, for each generator p/~ with u G 14 in the same 
side of I as pi (respectively pj), we consider the triplet (e, j, k) (respectively (e, i, k)). 
Again resorting to eq. H5.6J) and Lemma 13.81 we conclude that pk is also centered. 
Finally, if a generator pk with v G Vk is such that pk G I, any of the triplets (e, j, k) or 
(e, i, £;) can be invoked in a similar argument to ensure that pk is centered. Therefore, 
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P G A DC (Q), and hence (Z_ Ln (aig v ),w DC H Q n ) C A DC (Q). □ 

Proposition 5.8. For £/ie dynamical system l|5.5|l . i/ie generators' location P = 

(p±, . .. ,Pn) converges asymptotically to the largest weakly invariant set contained in 

the closure of A SP (Q) = {P G Q n \ i G I(P) => Pi G IC(V$)}. 

Proof. Let a G £Ln(9sm v )^fep(P)- By definition, a = Ln[smy](P) ■ Ci f° r ai l 

C £ dHsp(P)- Let e G Edsp(V(P)). If e is an edge of type (a), i.e. a segment of the 

bisector determined by pi and P j, we compute (cf. Proposition ^. 4J) . 



From Lemma I3.8r iii'> we conclude that a > 0, and the inequality is strict if either 
Pi G - IC(Vi) or pj £ IC(Vj). The same conclusion can be derived if e is a ver- 
tex of type (b). Therefore, max£L n (c) smv ) Hsp(P) > H SP (P) = 0. 
Now, resorting to the LaSalle principle ( Theorem 12. 9fl . we deduce that the solution 
P : [0, +oo) — > Q n starting from Pq converges to the largest weakly invariant set con- 
tained in ^Ln(9sm v ),Wsp 

nWgp(< H SP (P ) 1 Po)r\Q n . From cq. l|577|l. and resorting to 
Proposition 14.61 and Lemma 13.81 one can also show that ^Ln(9sm v ),'Hsp ^ Q n ^ s ec l ua l 
to A SP (Q). □ 

Remark 5.9. The sets Adc( ( 5) an d A$p(Q) are not closed in general. If dim Q = 
1, then it can be seen that they indeed arc. In higher dimensions one can find sequences 
{Pk G Q n | k G N} in these sets which converge to configurations P where not all 
active generators are centered. 

5.3. Distributed dynamical systems based on geometric centering. Here, 
we propose alternative distributed dynamical systems for the multi-center functions. 
Our design is directly inspired by the results in Theorems 14.71 and 14.81 on the critical 
points of the multi-center functions H.p>c an d Ti-sp- For i G {1, . . . , n}, consider the 
dynamical systems 



Alternatively, we may write P = CC(V(P)) — P and P G IC(V(P)) - P. Note that 
both systems are Voronoi-distributcd. Also, note that the vector field (|5.8J) is contin- 
uous, since the circumccnter of a polygon depends continuously on the location of its 
vertexes, and the location of the vertexes of the Voronoi partition depends continu- 
ously on the location of the generators; see [21] ■ However, eq. H5.9(l is a differential 
inclusion, since the incenter sets may not be singletons. By Lemma 12.71 the existence 
of solutions to eq. (|5.9J) is guaranteed by the following result. 

PROPOSITION 5.10. Consider the set-valued map IC(V) - Id : Q n — > 2^ K given 
by P i — ► IC(V(P)) — P. Then IC(V) — Id is upper semicontinuous with nonempty, 
compact and convex values. 

Proof. Clearly, the map IC(V) — Id takes nonempty and compact values. From 
Lemma l3.2l we also know that it takes convex values. Furthermore, since the identity 
map is continuous, it suffices to check that P i— > IC(V(P)) is upper semicontinuous. 
We then have to verify that, given Pq G Q n , for each e > 0, there exists 6 > such 
that 



a = Ln(<9smy)(P) • d e Fi 
= Ln(dsm Vi (p)){P) ■ n{d e Fi) + Ln(dsmy j( p))(P) ■ ir 3 (d e Fi) . 



(5.7) 



p i = CC(V i )-p i 
Pi elC(Vi)-Pi. 



(5.8) 
(5.9) 



IC(V(P)) ClC(V(P ))+P 2 „(0,e), if ||P - Poll < 5. 



(5.10) 
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Now, for each i, if lC(Vi(Po)) is not a singleton, then it is a segment (cf. Lemma Rj.2fl 
whose extremal points qn(Po), toC^o) are the intersection points of some bisectors 
of the edges of the Voronoi cell. It is clear that qi a (P) — > qi a (Po) when P — > Pq for 
a = 1,2. Therefore, given e > 0, one can choose > such that if \\P — Poll < 
<5j, then ||q iC( (P) — <7i Q (Po)|| < Since IC(V^(P)) is contained in the segment 

joining q a (P) and g i2 (P), we deduce IC(V Z (P)) C IC(^(P )) + P 2 (0,e/n). On the 
other hand, if lC(Vi (Po)) is a singleton, then it coincides with the intersection points 
qn{Po)i ■ ■ ■ ,Qim(Po) of some bisectors of the edges of the Voronoi cell. The above 
reasoning also guarantees that there exits Si > such that qi a (P) £ IC(V,;(Po)) + 
P 2 (0,e/n), a = l,...,m, if ||P — Po|| < Si. Since lC(Vi(P)) is contained in one 
of the segments joining the points qn(P), ■ ■ ■ , qim(P), we again deduce IC(Vi(P)) C 
IC(Vi(Po)) + B 2 (0, e/n). The statement in <|5.1U[) follows by taking the minimum of 
Si,...,S n . □ 

Having established the existence of solutions, one can also see that the compact 
set Q n is strongly invariant for the vector field CC(V) — Id and for the differential in- 
clusion IC(V) — Id. Next, we characterize the asymptotic convergence of the dynamical 
systems under study. 

Proposition 5.11. For the dynamical system Ij5.8|l (respectively (|5.9J) ). the gen- 
erators' location P = (pi, . . . ,p n ) converges asymptotically to the largest weakly invari- 
ant set contained in the closure of Adq(Q) (respectively in the closure of Asp(Q)). 

Proof. The proof of this result is parallel to the proofs of Propositions l5.7l and l5.8l 
The sequence of steps is the same as before, though now one resorts to Lemma rOT ii'l 
and Lemma I3.8f iv). The only additional observation is that, when computing the 
set- valued Lie derivative for eq. (|5.9|l . one has that a £ £ic(V)-id^Sp(-P) if an d only 
if there exists x £ IC(V(P)) such that a = (x - P) ■ (, for any ( £ dH S p(P). The 
application of Lemma |3. 81 guarantees that a > 0, and that the inequality is strict if 
any of the active generators is not in its corresponding incenter set. □ 

5.4. Simulations. To illustrate the performance of the distributed coordination 
algorithms, we include some simulation results. The algorithms are implemented in 
Mathematica as a single centralized program. We compute the bounded Voronoi dia- 
gram of a collection of points using the Mathematica package ComputationalGeometry 
We compute the circumcenter of a polygon via the algorithm in |28j and the incenter 
set via the LinearProgramming solver in Mathematica. Measuring displacements in 
meters, we consider the domain determined by the vertexes 

{(0, 0), (2.5, 0), (3.45, 1.5), (3.5, 1.6), (3.45, 1.7), (2.7, 2.1), (1., 2.4), (.2, 1.2)}. 

In Figs. 15.21 and 15.31 we illustrate the performance of the dynamical systems (|5.4H 
and (|5.8|l . respectively, minimizing the multi-circumcenter function "Hue- I n Figs. 15.41 
and 15.51 we illustrate the performance of the dynamical systems ()5.5|l and H5.9(l , 
respectively, maximizing the multi- incenter function Hsp ■ Observing the final config- 
urations in the four figures, one can verify, visually and numerically, that the active 
generators are asymptotically centered as forecast by our analysis. 

6. Conclusions. We have introduced two multi-center functions that provide 
quality-of-service measures for mobile networks. We have shown that both functions 
are globally Lipschitz and regular, and we have computed their generalized gradients. 
Furthermore, under certain technical conditions, we have characterized via nonsmooth 
analysis their critical points as center Voronoi configurations and as solutions of disk- 
covering and sphere-packing problems. We have also considered various algorithms 
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Fig. 5.2. "Toward the furthest" algorithm for 16 generators in a convex polygonal environment. 
The left (respectively, right) figure illustrates the initial (respectively, final) locations and Voronoi 
partition. The central figure illustrates the network evolution. After 2 seconds, the multi-center 
function is approximately .39504 meters. 




Fig. 5.3. "Move-toward-the-circumcenter" algorithm for 16 generators in a convex polygonal 
environment. The left (respectively, right) figure illustrates the initial (respectively, final) locations 
and Voronoi partition. The central figure illustrates the network evolution. After 20 seconds, the 
multi- center function is approximately 0.43273 meters. 

that extremize the multi-center functions. First, we considered the nonsmooth gra- 
dient flows induced by their respective generalized gradients. Second, we devised a 
novel strategy based on the generalized gradients of the 1-center functions of each 
generator. Third, we introduced and characterized a geometric centering strategy 
with resemblances to the classical Lloyd algorithm. We have unveiled the remarkable 
geometric interpretations of these algorithms, discussed their distributed character 
and analyzed their asymptotic behavior using nonsmooth stability analysis. 

Future directions of research include: (i) sharpening the asymptotic convergence 
results for the proposed dynamical systems, (ii) considering the setting of convex 
poly topes in R , for N > 2, (iii) understanding in what sense the proposed multi- 
circumcenter and the multi-incenter problems can be shown to be dual, and (iv) 
analyzing other meaningful geometric optimization problems and their relations with 
cooperative behaviors. 
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Symbol Description and page(s) when applicable 

Ap>c(Q) Set of configurations P £ Q n where all active generators are in 

the circumcenter of its own Voronoi region, 25 
Asp(Q) Set of configurations P £ Q n where all active generators are in 

the incenter set of its own Voronoi region, 26 
CC(Q) Circumcenter of polytope Q, 6 
CR(Q) Circumradius of polytope Q, 6 

Ds Distance function to the convex set S, 4 
Ed(Q) Edges of polygon Q, 4 
Edsp(V(P)) Edges where the value of TCsp(P) is attained, 6 

e(i) Edge of V(P) belonging to Vi and to the boundary of Q, 5 
e(i,j) Edge of V(P) determined by pi and pj, 5 
Fi (P) Smallest distance from pi to the boundary of V, (P) , 14 
Gi(P) Largest distance from pi to the boundary of Vi(P), 14 
df Generalized gradient of the locally Lipschitz function /, 7 
Tiuc Multi-circumcentcr function, 6 
Tisp Multi-inccnter function, 6 
IC(Q) Incenter set of polytope Q, 6 
IR(Q) Inradius of polytope Q, 6 

K[X] Filippov mapping associated with a measurable and essentially 
locally bounded mapping X : R N -> R N , 8 
X(e,i,j) Scalar function associated with the vertex v(e,i,j), 16 
Ln(5) Least-norm element of the convex set S, 7 
lgg(p) Largest distance from p to the boundary of Q, 10 
n{i,j, k) Scalar function associated with the vertex v(i,j, k), 17 
Af(P, i), Af(i) Set of neighbors of the ith. generator at configuration P, 4 
n e(i,j) Unit normal to e(i,j) pointing toward mt{Vi(P)) , 5 
H e (j) Unit normal to e(i) pointing toward int(Q), 5 
proj s Orthogonal projection onto the convex set S, 4 

TTi Canonical projection from Q n onto the ith factor, 4 
Cxf Set-valued Lie derivative of / with respect to X, 8 
smg(p) Smallest distance from p to the boundary of Q, 10 
v(i,j, k) Vertex of V(P) determined by pi, pj and pk, 4 
v(e, Vertex of V(P) determined by e € Ed(Q) and pi, pj, 4 
v(e, /, i) Vertex of V(P) determined by e, / £ Ed(Q) and p^ 4 
VeDc(V(P)) Vertcxes of V(P) where the value of Hdc(P) is attained, 6 
vrs(i;) Unit vector in the direction of ^ v £ R w , 4 
Vc(Q) Vertcxes of polygon Q, 4 
V(P) Voronoi partition of Q generated by P = (pi, . . . ,p n ), 4 
Zx.f Set formed by points x £ M. N such that belongs to Cxf{x), 9 



